#!/usr/bin/perl

# Script: wvr_ops.pl
#
# This script does everything necessary to construct the WVR phase corrections
# for a particular observation. It only requires the set of RPFITS files that
# need correcting, or the Miriad dataset made from those RPFITS files,
# and it will do all the rest automatically.
#
# Copyright 2013 Jamie Stevens
# Email: Jamie.Stevens@csiro.au
#
# This program is free software: you can redistribute it and/or modify
# it under the terms of the GNU General Public License as published by
# the Free Software Foundation, either version 3 of the License, or
# (at your option) any later version.
#
# This program is distributed in the hope that it will be useful,
# but WITHOUT ANY WARRANTY; without even the implied warranty of
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
# GNU General Public License for more details.
#
# You should have received a copy of the GNU General Public License
# along with this program.  If not, see <http://www.gnu.org/licenses/>.
#

# EDIT: POINT THIS TO THE LOCATION OF YOUR ATNF::MoniCA LIBRARY.
# IF YOU HAVE IT INSTALLED SYSTEM-WIDE, YOU CAN COMMENT OUT THIS LINE.
use lib "/n/ste616/src/open-monica/clientlibs/perl";

use Astro::Time;
use ATNF::MoniCA;
use PDL;
use PDL::GSL::INTERP;
use PDL::Graphics::PGPLOT;
use PGPLOT;
use strict;

# The arguments should be a list of RPFITS files that need corrections.

print STDERR "START: wvr_ops.pl\n";

my @rpfits;
my @fpaths;
my $data_is_miriad = 1;
my ($mintime, $maxtime);
my $refresh = 0;
my $cycle_time = 5;
my $align_wvr = 0;
my $verbose = 0;
my @excluded_ants;
my $blcutoff;
for (my $i=0; $i<=$#ARGV; $i++) {
    # This should be an RPFITS file, check that it exists.
    if ($ARGV[$i] eq "--cycle" ||
	$ARGV[$i] eq "-c") {
	# The user has set the cycle time.
	$i++;
	$cycle_time = $ARGV[$i];
    } elsif ($ARGV[$i] eq "--refresh") {
	# Force reloading of all information.
	$refresh = 1;
    } elsif ($ARGV[$i] eq "--wvr-align") {
	# Align the times to the WVR data.
	$align_wvr = 1;
    } elsif ($ARGV[$i] eq "--verbose" ||
	     $ARGV[$i] eq "-v") {
	# Be more verbose.
	$verbose = 1;
    } elsif ($ARGV[$i] eq "--rpfits") {
	# Use RPFITS files instead of a Miriad dataset.
	$data_is_miriad = 0;
#    } elsif ($ARGV[$i] eq "--cutoff") {
#	# Zero correction on baselines shorter than this specified length (m).
#	$i++;
#	$blcutoff = $ARGV[$i];
    } elsif ($ARGV[$i] eq "--exclude" ||
	     $ARGV[$i] eq "-x") {
	# Exclude an antenna from the corrections.
	$i++;
	push @excluded_ants, $ARGV[$i];
    } elsif ($ARGV[$i] eq "--help" ||
	     $ARGV[$i] eq "-h") {
	# Print out the help message and exit.
	&print_help();
    } elsif (-e $ARGV[$i]) {
	push @rpfits, $ARGV[$i];
    }
}

# The available baselines.
my @baselines = ( [2, 1], [3, 1], [3, 2], [4, 1], [4, 2],
		  [4, 3], [5, 1], [5, 2], [5, 3], [5, 4],
		  [6, 1], [6, 2], [6, 3], [6, 4], [6, 5] );
my @baseline_lengths;
#my @baseline_exclude;
my @mons = qw/JAN FEB MAR APR MAY JUN JUL AUG SEP OCT NOV DEC/;


# Get the date range to produce the WVR corrections for in day fractions.
# Get some information from our RPFITS files or Miriad datasets.
# Time range and frequencies using the atcacode rpfits_summary or
# Miriad's uvlist.
if ($#rpfits<0) {
    my $type = "RPFITS files";
    if ($data_is_miriad) {
	$type = "Miriad datasets";
    }
    print STDERR "FATAL: No ".$type." were given!\n";
    print STDERR "       Please specify the ".$type." containing the ".
	"data that you wish to correct.\n";
    exit;
}

my %freqlist;
my %sourcelist;
my @obsorder;

my ($rpfread, $uvindex);
if (!$data_is_miriad) {
    $rpfread = "rpfits_summary";
    if (&check_path($rpfread) == 0) {
	die "FATAL: Require atcacode's rpfits_summary to be available and in the PATH.\n";
    }
    print STDERR "STATE: Getting information from RPFITS files...\n";
} else {
    $uvindex = "uvindex";
    if (&check_path($uvindex) == 0) {
	die "FATAL: Require Miriad's uvsplit to be available and in the PATH.\n";
    }
    print STDERR "STATE: Getting information from Miriad datasets...\n";
}

for (my $i=0; $i<=$#rpfits; $i++) {
    if ($verbose) {
	if (!$data_is_miriad) {
	    print STDERR "       Reading file ".$rpfits[$i]."... ";
	} else {
	    print STDERR "       Reading dataset ".$rpfits[$i]."... ";
	}
    }
    my ($year, $month, $date, $ut, $fpath, $command);
    $fpath = "";
    if (!$data_is_miriad) {
	if ($rpfits[$i] =~ /^(.*\/*)(\d{4})\-(\d{2})\-(\d{2})\_(\d{4})\.(.*)$/) {
	    # This looks correct.
	    $year = $2;
	    $month = $3;
	    $date = $4;
	    $fpath = $1;
	    $rpfits[$i] = $2."-".$3."-".$4."_".$5.".".$6;
	} else {
	    if ($verbose) {
		print STDERR "Not an RPFITS file!\n";
	    }
	    splice @rpfits, $i, 1;
	    $i--;
	    next;
	}
	$command = $rpfread." ".$fpath.$rpfits[$i];
	$fpaths[$i] = $fpath;
    } else {
	# Remove any trailing /.
	$rpfits[$i] =~ s/\/$//;
	my $fpath = "";
	if ($rpfits[$i] =~ /^(.*)\/(.*)$/) {
	    # There's a path here.
	    $fpath =~ $1."/";
	    $rpfits[$i] = $2;
	}
	if (!-d $fpath.$rpfits[$i] ||
	    !-e $fpath.$rpfits[$i]."/visdata") {
	    if ($verbose) {
		print STDERR "Not a Miriad dataset!\n";
	    }
	    splice @rpfits, $i, 1;
	    $i--;
	    next;
	}
	$command = $uvindex." vis=".$fpath.$rpfits[$i];
	$fpaths[$i] = $fpath;
    }
    my $cache_file = $fpath.".cache.".$rpfits[$i];
    my $from_cache = 0;
    if (-e $cache_file && $refresh == 0) {
	open(RPF, $cache_file);
	$from_cache = 1;
    } else {
	open(RPF, "-|") || exec $command;
	open(CACHE, ">".$cache_file);
    }
    my $oldut = 0;
    my $fdo = 0;
    while(<RPF>) {
	chomp;
	my $line = $_;
	my @els = split(/\s+/, $line);
	if ($from_cache == 0) {
	    print CACHE $line."\n";
	}
	if (!$data_is_miriad && $line =~ /^\>\>\>/) {
	    # We don't care about this line.
	    next;
	}
	if (!$data_is_miriad ||
	    $els[0] =~ /^(\d{2})(\D{3})(\d{2})\:(\d{2})\:(\d{2})\:(.+)$/) {
	    if (!$data_is_miriad) {
		# Get rid of all the parentheses.
		$line =~ s/[\(\)]//g;
		@els = split(/\s+/, $line);
		($year, $month, $date) = split(/\-/, $els[1]);
		$ut = $els[2] / 86400.0;
	    } else {
		($year, $month, $date) = ($1 + 2000, str2month($2), $3);
		my ($h, $m, $s) = ($4, $5, $6);
		$ut = hms2time($h, $m, $s);
	    }
	    my $mjd = cal2mjd($date, $month, $year, $ut);
	    if ($mintime == 0) {
		$mintime = $mjd;
		$maxtime = $mjd;
	    } else {
		$mintime = ($mjd < $mintime) ? $mjd : $mintime;
		$maxtime = ($mjd > $maxtime) ? $mjd : $maxtime;
	    }
	    # Get the source information.
	    my $src;
	    if (!$data_is_miriad) {
		$src = $els[3];
	    } else {
		$src = $els[1];
	    }
	    if (!defined $sourcelist{$src}) {
		$sourcelist{$src} = {
		    'start_times' => [ $mjd ]
		};
	    } else {
		push @{$sourcelist{$src}->{'start_times'}}, $mjd;
	    }
	    push @obsorder, {
		'source' => $src,
		'start_time' => $mjd
	    };
	}

	if (!$data_is_miriad || $fdo) {
	    if ($data_is_miriad && $line eq "") {
		$fdo = 0;
		next;
	    }
	    if ($data_is_miriad && $fdo == 1) {
		$fdo++;
		next;
	    }
	    my ($tfreq, $bandw, $cfreq);
	    if (!$data_is_miriad) {
		# Go through the frequencies listed.
		for (my $f=4; $f<=$#els; $f+=3) {
		    $tfreq = "f".$els[$f];
		    $bandw = $els[$f + 1];
		    $cfreq = $els[$f];
		    if (!defined $freqlist{$tfreq}) {
			$freqlist{$tfreq} = {
			    'centre' => $cfreq,
			    'bandwidths' => [ $bandw ]
			};
		    } else {
			&unique_push($freqlist{$tfreq}->{'bandwidths'}, $bandw);
		    }
		}
	    } else {
		my $tchan = ($els[1] % 2) ? ($els[1] - 1) : $els[1];
		my $tsfreq = $els[2];
		my $tschan = $els[3];
		$cfreq = (($tchan / 2) * $tschan + $tsfreq) * 1000.0;
		$bandw = $tchan * $tschan * 1000.0;
		$tfreq = sprintf "f%.0f", $cfreq;
		if (!defined $freqlist{$tfreq}) {
		    $freqlist{$tfreq} = {
			'centre' => $cfreq,
			'bandwidths' => [ $bandw ]
		    };
		} else {
		    &unique_push($freqlist{$tfreq}->{'bandwidths'}, $bandw);
		}
	    }
	}
	if ($data_is_miriad && $line =~ /^Frequency Configuration \d+/) {
	    $fdo = 1;
	}
    }
    close(RPF);
    if ($from_cache == 0) {
	close(CACHE);
    }
    if ($verbose) {
	print STDERR "Done.\n";
    }
}
if (!$data_is_miriad) {
    print STDERR "STATE: Finished reading RPFITS files.\n";
} else {
    print STDERR "STATE: Finished reading Miriad datasets.\n";
}

if ($#rpfits<0) {
    my $type = "RPFITS files";
    if ($data_is_miriad) {
	$type = "Miriad datasets";
    }
    print STDERR "FATAL: No valid ".$type." were given!\n";
    print STDERR "       Please specify the ".$type." containing the ".
	"data that you wish to correct.\n";
    exit;
}

my @lofreqs;
my @hifreqs;
my @cefreqs;
my @cewaves;
my $maxwave = 0;
foreach my $flabel (keys %freqlist) {
    for (my $i=0; $i<=$#{$freqlist{$flabel}->{'bandwidths'}}; $i++) {
	my $l = $freqlist{$flabel}->{'centre'} -
	    ($freqlist{$flabel}->{'bandwidths'}->[$i] / 2);
	my $h = $freqlist{$flabel}->{'centre'} +
	    ($freqlist{$flabel}->{'bandwidths'}->[$i] / 2);
	if ($verbose) {
	    print STDERR "INFO: found freq range $l -- $h (MHz).\n";
	}
	push @lofreqs, $l;
	push @hifreqs, $h;
	push @cefreqs, $freqlist{$flabel}->{'centre'};
	my $lambda = 299792458000 / ($freqlist{$flabel}->{'centre'} * 1E6);
	if ($lambda > $maxwave) {
	    $maxwave = $lambda;
	}
	push @cewaves, $lambda;
    }
}

# Fill in the baseline lengths.
print STDERR "STATE: Determining baseline lengths...";
my @axpos;
my @aypos;
my @azpos;
if ($data_is_miriad) {
    # Get the antenna positions from Miriad.
    my $bcmd = "uvlist vis=".$fpaths[0].$rpfits[0].
	" options=array";
    open(BASE, "-|") || exec $bcmd;
    while(<BASE>) {
	chomp;
	my @els = split(/\s+/);
	if ($els[1] =~ /^\d$/) {
	    $axpos[$els[1] - 1] = $els[2];
	    $aypos[$els[1] - 1] = $els[3];
	    $azpos[$els[1] - 1] = $els[4];
	}
    }
    close(BASE);
} else {
    # Get the antenna positions from the RPFITS header.
    my $bcmd = "head -n 1 ".$fpaths[0].$rpfits[0]." | ".
	"strings | fold | grep CA0";
    open(BASE, "-|") || exec $bcmd;
    while(<BASE>) {
	chomp;
	my @els = split(/\s+/);
	if ($els[1] =~ /^\d$/) {
	    $axpos[$els[1] - 1] = $els[4];
	    $aypos[$els[1] - 1] = $els[5];
	    $azpos[$els[1] - 1] = $els[6];
	}
    }
    close(BASE);
}
for (my $i=0; $i<=$#baselines; $i++) {
    $baseline_lengths[$i] = sqrt(
	($axpos[$baselines[$i]->[0] - 1] - $axpos[$baselines[$i]->[1] - 1])**2 +
	($aypos[$baselines[$i]->[0] - 1] - $aypos[$baselines[$i]->[1] - 1])**2 +
	($azpos[$baselines[$i]->[0] - 1] - $azpos[$baselines[$i]->[1] - 1])**2);
#    if ($baseline_lengths[$i] < $blcutoff) {
#	push @baseline_exclude, $i;
#    }
    if ($verbose) {
	printf STDERR "\n       CA0%1d-CA0%1d = %.3f m",
	$baselines[$i]->[0], $baselines[$i]->[1], $baseline_lengths[$i];
#	if ($baseline_exclude[$#baseline_exclude] == $i) {
#	    printf STDERR " (excluded)";
#	}
    }
}
if ($verbose) {
    print STDERR "\nSTATE: Baseline lengths determined.\n";
} else {
    print STDERR " done.\n";
}

# Now get the MoniCA data for the WVR for this time range.
# The MoniCA points:
my $allpointsrequired = 0;
my @systems = ( "SkyRF1", "SkyRF2", "SkyRF3", "SkyRF4" );
my @systems_check = ( "PaddleRF1", "PaddleRF2", "PaddleRF3", "PaddleRF4" );
if ($allpointsrequired == 1) {
    push @systems, ( "SkyRF5", "T1", "T1a", "T2", "T3", "T4",
		     "T5", "Irf", "Ish", "+10VRef", "-10VRef" );
}
my $point_prefix = "receiver.WVR";

# Connect to the MoniCA server.
my $monserver = "monhost-nar.atnf.csiro.au";
print STDERR "STATE: Connecting to MoniCA server ".$monserver."...";

my $monhandle = monconnect($monserver);
if (!defined $monhandle) {
    print STDERR "\n";
    print STDERR "FATAL: Unable to connect to the MoniCA server.\n";
    print STDERR "       Ensure that you are able to access the CASS ".
	"internal network from this machine.\n";
#    print STDERR " failed!\n";
    exit;
}
print STDERR " done.\n";

# Get the dUTC.
if ($verbose) {
    print STDERR "STATE: Getting dUTC...";
}
my $dUTC = monpoll($monhandle, "caclock.misc.clock.dUTC");
if ($verbose) {
    print STDERR " ".$dUTC->val()." secs.\n";
}

# Get all the data.
# Search the time range plus/minus one hour.
my $minsearch = $mintime - 1/24;
my $maxsearch = $maxtime + 1/24;
my %wvrdata;
print STDERR "STATE: Getting WVR points between ".
    mjd2time($minsearch)." -- ".mjd2time($maxsearch)."...";
if ($verbose) {
    print STDERR "\n";
}
for (my $p=0; $p<=$#systems; $p++) {
    my $psuffix = $point_prefix.".".$systems[$p];
    if ($verbose) {
	print STDERR "       ca0?.".$psuffix."...";
    }
    $wvrdata{$psuffix} = [];
    for (my $a=1; $a<=6; $a++) {
	# Loop over antenna.
	if ($verbose) {
	    print STDERR " ca0".$a;
	}
	my $tpoint = "ca0".$a.".".$psuffix;
	$wvrdata{$psuffix}->[$a - 1] = [];
	my $done = 0;
	my @data;
	my $cache_file = ".cache.".$psuffix.".ca0".$a;
	if (-e $cache_file && $refresh == 0) {
	    open(CACHE, $cache_file);
	    while(<CACHE>) {
		chomp;
		push @data, new MonBetweenPoint($_);
	    }
	    close(CACHE);
	} else {
	    @data = monbetween($monhandle, $minsearch, $maxsearch, $tpoint);
	    open(CACHE, ">".$cache_file);
	    for (my $i=0; $i<=$#data; $i++) {
		print CACHE $data[$i]->bat."\t".$data[$i]->val."\t".
		    $data[$i]->errorstate."\n";
	    }
	    close(CACHE);
	}	    
	push @{$wvrdata{$psuffix}->[$a - 1]}, @data;
	if ($verbose) {
	    print STDERR "(".($#data + 1).")";
	}
    }
    if ($verbose) {
	print STDERR " done.\n";
    }
}
if ($verbose) {
    print STDERR "STATE: All points obtained.\n";
} else {
    print STDERR " done.\n";
}
print STDERR "STATE: Checking WVR state for this observation... ";
if ($verbose) {
    print STDERR "\n";
}
for (my $p=0; $p<=$#systems_check; $p++) {
    my $psuffix = $point_prefix.".".$systems_check[$p];
    if ($verbose) {
	print STDERR "       ca0?.".$psuffix."...";
    }
    $wvrdata{$psuffix} = [];
    for (my $a=1; $a<=6; $a++) {
	# Loop over antenna.
	if ($verbose) {
	    print STDERR " ca0".$a;
	}
	my $tpoint = "ca0".$a.".".$psuffix;
	$wvrdata{$psuffix}->[$a - 1] = [];
	my $done = 0;
	my @data;
	my $cache_file = ".cache.".$psuffix.".ca0".$a;
	if (-e $cache_file && $refresh == 0) {
	    open(CACHE, $cache_file);
	    while(<CACHE>) {
		chomp;
		push @data, new MonBetweenPoint($_);
	    }
	    close(CACHE);
	} else {
	    @data = monbetween($monhandle, $minsearch, $maxsearch, $tpoint);
	    open(CACHE, ">".$cache_file);
	    for (my $i=0; $i<=$#data; $i++) {
		print CACHE $data[$i]->bat."\t".$data[$i]->val."\t".
		    $data[$i]->errorstate."\n";
	    }
	    close(CACHE);
	}	    
	push @{$wvrdata{$psuffix}->[$a - 1]}, @data;
	if ($verbose) {
	    print STDERR "(".($#data + 1).")";
	}
    }
    if ($verbose) {
	print STDERR " done.\n";
    }
}
if ($verbose) {
    print STDERR "STATE: All points obtained.\n";
} else {
    print STDERR " done.\n";
}
# Output some warnings if required.
my @nopaddles;
for (my $a=1; $a<=6; $a++) {
    my $total_paddles = 0;
    for (my $i=0; $i<=$#systems_check; $i++) {
	my $csuffix = $point_prefix.".".$systems_check[$i];
	$total_paddles += ($#{$wvrdata{$csuffix}->[$a - 1]} + 1);
    }
    if ($total_paddles == 0) {
	push @nopaddles, $a;
    }
}
if ($#nopaddles >= 0) {
    print STDERR "WARNING: WVR Paddles were not functioning on";
    for (my $i=0; $i<=$#nopaddles; $i++) {
	print STDERR " CA0".$nopaddles[$i];
    }
    print STDERR " during these observations.\n";
    print STDERR "WARNING: Corrections may not be accurate.\n";
}

# Make our piddles for later.
my @wvrpiddles;
my ($common_min_day, $common_max_day);
print STDERR "STATE: Converting to piddles...";
my $sa = 0;
for (my $a=1; $a<=6; $a++) {
    my $vref = $wvrdata{$point_prefix.".".$systems[0]}->[$a - 1];
    my $tmp = {
	'Day' => []
    };
    for (my $p=0; $p<=$#systems; $p++) {
	$tmp->{$systems[$p]} = [];
    }
    for (my $v=0; $v<=$#{$vref}; $v++) {
	# Make the MJD array.
	push @{$tmp->{'Day'}}, bat2mjd($vref->[$v]->bat);
	# And the arrays for all the other systems.
	for (my $p=0; $p<=$#systems; $p++) {
	    my $tref = $wvrdata{$point_prefix.".".$systems[$p]}->[$a - 1];
	    push @{$tmp->{$systems[$p]}}, $tref->[$v]->val;
	}
    }
    # Turn each array now into a piddle.
    $wvrpiddles[$a - 1] = {};
    foreach my $s (keys %{$tmp}) {
	$wvrpiddles[$a - 1]->{$s} = pdl $tmp->{$s};
    }
    # Get the minimum and maximum times returned from MoniCA.
    my ($mn, $mx) = minmax($wvrpiddles[$a - 1]->{'Day'});
    if ($sa == 0 && $mn > 0 && $mx > 0) {
	$common_min_day = $mn;
	$common_max_day = $mx;
	$sa = 1;
    } else {
	# We need to find the largest time range that all antenna
	# have in common, that's why this next block appears backwards.
	$common_min_day = ($mn > $common_min_day) ? $mn : $common_min_day;
	$common_max_day = ($mx < $common_max_day) ? $mx : $common_max_day;
    }
}
print STDERR " done.\n";

# Now adjust the minimum and maximum times to lie within the available
# time range, but remain on the cycle time boundaries defined by the
# RPFITS files.
my $cycletime_indays = $cycle_time / 86400;
if ($align_wvr) {
    $mintime = $common_min_day;
    $maxtime = $common_max_day;
} else {
    while ($mintime > ($common_min_day + $cycletime_indays)) {
	$mintime -= $cycletime_indays;
    }
    while ($maxtime < ($common_max_day - $cycletime_indays)) {
	$maxtime += $cycletime_indays;
    }
}
$mintime = POSIX::ceil($mintime * 86400) / 86400 + $cycletime_indays;
$maxtime = POSIX::ceil($maxtime * 86400) / 86400 - $cycletime_indays;

# Make some files for legacy applications.
print STDERR "STATE: Creating legacy files...";
for (my $a=1; $a<=6; $a++) {
    my $file = "ca0".$a;
    open(OUT, ">".$file);
    foreach (0 .. $wvrpiddles[$a - 1]->{'Day'}->nelem - 1) {
	print OUT $wvrpiddles[$a - 1]->{'Day'}->index($_);
	for (my $p=0; $p<=$#systems; $p++) {
	    print OUT " ".$wvrpiddles[$a - 1]->{$systems[$p]}->index($_);
	}
	print OUT "\n";
    }
    close(OUT);
}
print STDERR " done.\n";

# Some constants.
my $rhowater = 1e6;
my $D = 76000.0 / 217.0 / 0.001763 / 100.0; # Scale height.

# Standard atmosphere vapour paths.
my $dTdL16 = 0.043209199;
my $dTdL18 = 0.089893485;
my $dTdL22 = 0.234429107;
my $dTdL25 = 0.16258403;

my $Weight16 = 0.020440843;
my $Weight18 = 0.088471382;
my $Weight22 = 0.601685362;
my $Weight25 = 0.289402413;

# Calibration factors for each WVR.
my @calibration_data = (
    # Unit 1, CA01
    [ { 'a' => 77.1384, 'b' => -319.2057 },
      { 'a' => 76.9287, 'b' => -251.4061 },
      { 'a' => 73.6508, 'b' => -249.9106 },
      { 'a' => 95.7108, 'b' => -316.0245 } ],
    # Unit 2, CA02
    [ { 'a' => 147.9768, 'b' => -370.3523 },
      { 'a' => 72.4621, 'b' => -266.8057 },
      { 'a' => 70.0063, 'b' => -266.8866 },
      { 'a' => 127.5982, 'b' => -334.6185 } ],
    # Unit 3, CA03
    [ { 'a' => 97.6599, 'b' => -300.8925 },
      { 'a' => 167.5655, 'b' => -267.0325 },
      { 'a' => 70.8773, 'b' => -287.0931 },
      { 'a' => 117.257, 'b' => -329.4508 } ],
    # Unit 7, CA04
    [ { 'a' => 75.7153, 'b' => -305.0161 },
      { 'a' => 73.6768, 'b' => -251.8156 },
      { 'a' => 121.0595, 'b' => -320.542 },
      { 'a' => 133.8446, 'b' => -382.5208 } ],
    # Unit 5, CA05
    [ { 'a' => 67.3326, 'b' => -288.8005 },
      { 'a' => 71.3074, 'b' => -279.3064 },
      { 'a' => 76.5698, 'b' => -316.7808 },
      { 'a' => 98.4227, 'b' => -352.5096 } ],
    # Unit 6, CA06
    [ { 'a' => 72.8825, 'b' => -296.8122 },
      { 'a' => 70.2582, 'b' => -244.1924 },
      { 'a' => 101.1862, 'b' => -265.6255 },
      { 'a' => 117.1777, 'b' => -295.027 } ]
    );

# Form a master timeline from the first to the last MJD in cycle time increments.
my $num_points = ($maxtime - $mintime) * 86400 / $cycle_time;
my $master_timeline = sequence($num_points) * $cycle_time / 86400 + $mintime;
# CHECKED HERE, WORKS AS EXPECTED, MATCHES wvr_phase.pro
#print $master_timeline->index($_)."   " foreach (0 .. $master_timeline->nelem-1);
#print "\n";
# Set up all the required interpolations.
my @interps;
print STDERR "STATE: Setting up interpolaters...";
for (my $a=1; $a<=6; $a++) {
    print STDERR " ca0".$a;
    $interps[$a - 1] = {};
    $interps[$a - 1]->{'DayFrac'} = 
	PDL::GSL::INTERP->init('cspline', 
			       $wvrpiddles[$a - 1]->{'Day'},
			       $wvrpiddles[$a - 1]->{'Day'});
    for (my $p=0; $p<=$#systems; $p++) {
	$interps[$a - 1]->{$systems[$p]} = 
	    PDL::GSL::INTERP->init('cspline', $wvrpiddles[$a - 1]->{'Day'},
				   $wvrpiddles[$a - 1]->{$systems[$p]});
    }
}
print STDERR " done.\n";
# And now do the interpolation.
my @DayFrac;
my @C165;
my @C189;
my @C229;
my @C255;
print STDERR "STATE: Performing interpolations...";
for (my $a=1; $a<=6; $a++) {
    print STDERR " ca0".$a;
    my @tmp = ( [], [], [], [], [] );
    foreach (0 .. $master_timeline->nelem - 1) {
	my $ttime = $master_timeline->index($_);
	push @{$tmp[0]}, $interps[$a - 1]->{'SkyRF1'}->eval($ttime);
	push @{$tmp[1]}, $interps[$a - 1]->{'SkyRF2'}->eval($ttime);
	push @{$tmp[2]}, $interps[$a - 1]->{'SkyRF3'}->eval($ttime);
	push @{$tmp[3]}, $interps[$a - 1]->{'SkyRF4'}->eval($ttime);
	push @{$tmp[4]}, $interps[$a - 1]->{'DayFrac'}->eval($ttime);
    }
    $C165[$a - 1] = pdl $tmp[0];
    $C189[$a - 1] = pdl $tmp[1];
    $C229[$a - 1] = pdl $tmp[2];
    $C255[$a - 1] = pdl $tmp[3];
    $DayFrac[$a - 1] = pdl $tmp[4];
}
print STDERR " done.\n";
# CHECKED HERE, WORKS AS EXPECTED, MATCHES wvr_phase.pro
#print $C255[0]->index($_)."  " foreach (0 .. $C255[0]->nelem-1);
#print "\n";

# Convert the voltages to temperatures.
my @T165;
my @T189;
my @T229;
my @T255;
print STDERR "STATE: Converting to temperatures...";
for (my $a=1; $a<=6; $a++) {
    print STDERR " ca0".$a;
    $T165[$a - 1] = $calibration_data[$a - 1][0]->{'a'} * $C165[$a - 1] -
	$calibration_data[$a - 1][0]->{'b'};
    $T189[$a - 1] = $calibration_data[$a - 1][1]->{'a'} * $C189[$a - 1] -
	$calibration_data[$a - 1][1]->{'b'};
    $T229[$a - 1] = $calibration_data[$a - 1][2]->{'a'} * $C229[$a - 1] -
	$calibration_data[$a - 1][2]->{'b'};
    $T255[$a - 1] = $calibration_data[$a - 1][3]->{'a'} * $C255[$a - 1] -
	$calibration_data[$a - 1][3]->{'b'};
}
print STDERR " done.\n";
# CHECKED HERE, WORKS AS EXPECTED, MATCHES wvr_phase.pro
#print $T255[0]->index($_)." " foreach (0 .. $T255[0]->nelem-1);
#print "\n";

my $ua = 1;
while (&ant_excluded($ua)) {
    $ua++;
}
if ($verbose) {
    print "INFO: Using CA0".$ua." as \"reference\" antenna.\n";
}
my $df_offset = floor($obsorder[0]->{'start_time'});
for (my $fi=0; $fi<=$#cefreqs; $fi++) {
    my $devname = sprintf "WVR_Phase_%ds_%d_all.ps/cps", 
    $cycle_time, floor($cefreqs[$fi]);#, $wvrTime->index(0);
    my $offsets_plt = PDL::Graphics::PGPLOT::Window->new(
	'Device' => $devname,
	'NXPanel' => 1,
	'NYPanel' => 3);
    print STDERR "STATE: Calculating the WVR phase corrections at ".
	$cefreqs[$fi]." MHz ... ";
    open(OUT, ">atwvr_input_".$cycle_time."s.".floor($cefreqs[$fi]));
    #for (my $oo=11; $oo<12; $oo++) { # TEST MODE
    for (my $oo=0; $oo<$#obsorder; $oo++) {
	my $mjd_start = $obsorder[$oo]->{'start_time'};
	while ($obsorder[$oo + 1]->{'source'} eq
	       $obsorder[$oo]->{'source'}) {
	    $oo++;
	}
	my $mjd_end = $obsorder[$oo + 1]->{'start_time'};

	# Get the date range to produce the WVR corrections for in day fractions.
	my $df_start = $mjd_start - floor($mjd_start);
	my $df_end = $df_start + ($mjd_end - $mjd_start);
#    print "Start: ".$mjd_start." end: ".$mjd_end."\n";
#    print "Day fraction start: ".$df_start." end: ".$df_end."\n";

	# CHECKED HERE, MATCHES wvr_phase.pro

	# Now select the time frame we're interested in and store it.
	my $sel_data = $DayFrac[$ua - 1]->where($DayFrac[$ua - 1] >= $mjd_start &
						$DayFrac[$ua - 1] <= $mjd_end);
#    print $sel_data."\n";
	my $count = $sel_data->nelem;
#    print "Number of matching timestamps: ".$count."\n";
    
	if ($count == 0) {
#	print "No matching timestamps were found. Examine the data files ".
#	    "and make sure the MJD timestamps are within the range of ".
#	    "you observing data!\n";
#	exit;
	    next;
	}
    
	my $wvrTime;
	my @TT165;
	my @TT189;
	my @TT229;
	my @TT255;
	for (my $a=1; $a<=6; $a++) {
	    ($wvrTime, $TT165[$a - 1], $TT189[$a - 1],
	     $TT229[$a - 1], $TT255[$a - 1]) =
		 where($DayFrac[$ua - 1], $T165[$a - 1], $T189[$a - 1],
		       $T229[$a - 1], $T255[$a - 1],
		       $DayFrac[$ua - 1] >= $mjd_start &
		       $DayFrac[$ua - 1] <= $mjd_end);
	    $wvrTime = $wvrTime - $df_offset;
	}
	
#    my $xplt = PDL::Graphics::PGPLOT::Window->new(
#	'Device' => 'test1.ps/cps',
#	'NXPanel' => 2, 'NYPanel' => 3);
#    for (my $a=1; $a<=6; $a++) {
#	my $data = $TT229[$a - 1] - $TT229[5];
#	$xplt->line($wvrTime, $data);
#	$xplt->label_axes("Time since scan start [d]",
#			  "[V]", "Tsky at 22.9 GHz, BL= 6-".$a);
#    }
#    $xplt->close();
    # CHECKED HERE, VERY SIMILAR TO wvr_phase.pro

	# Calculate the WVR phase for every baseline.
	my @dT165;
	my @dT189;
	my @dT229;
	my @dT255;
	my @Wsum;
	my @WVRPhi;
	my @WVRPhi00;
	for (my $bl=0; $bl<=$#baselines; $bl++) {
	    # Find which antennas are involved in this baseline.
	    my $a1 = $baselines[$bl]->[0];
	    my $a2 = $baselines[$bl]->[1];
	    
	    # Begin calculating the DC offsets based on average
	    # subtraction over a given time period.
	    # The first antenna.
	    $TT165[$a1 - 1] -= median($TT165[$a1 - 1]);
	    $TT189[$a1 - 1] -= median($TT189[$a1 - 1]);
	    $TT229[$a1 - 1] -= median($TT229[$a1 - 1]);
	    $TT255[$a1 - 1] -= median($TT255[$a1 - 1]);
	    # The second antenna.
	    $TT165[$a2 - 1] -= median($TT165[$a2 - 1]);
	    $TT189[$a2 - 1] -= median($TT189[$a2 - 1]);
	    $TT229[$a2 - 1] -= median($TT229[$a2 - 1]);
	    $TT255[$a2 - 1] -= median($TT255[$a2 - 1]);
	    
	    # Build the difference between the units, and the
	    # vapour path. The dT arrays will not be a delta
	    # temperature.
	    $dT165[$bl] = ($TT165[$a2 - 1] - $TT165[$a1 - 1]) / $dTdL16;
	    $dT189[$bl] = ($TT189[$a2 - 1] - $TT189[$a1 - 1]) / $dTdL18;
	    $dT229[$bl] = ($TT229[$a2 - 1] - $TT229[$a1 - 1]) / $dTdL22;
	    $dT255[$bl] = ($TT255[$a2 - 1] - $TT255[$a1 - 1]) / $dTdL25;
	    
	    # And build the weighted sum.
	    $Wsum[$bl] = $dT165[$bl] * $Weight16 +
		$dT189[$bl] * $Weight18 +
		$dT229[$bl] * $Weight22 +
		$dT255[$bl] * $Weight25;
	    $WVRPhi[$bl] = 2 * 180.0 * $Wsum[$bl] / $cewaves[$fi];
	}

	# Set those baselines to zero that are shorter than
	# the cutoff length.
#	for (my $ce=0; $ce<=$#baseline_exclude; $ce++) {
#	    $WVRPhi[$baseline_exclude[$ce]] *= 0.0;
#	}

	my $plt_x = $wvrTime * 1440.0;
	my @plot_groups = (
	    [ 10, 11, 12, 13, 14 ],
	    [ 0, 1, 2, 3, 4, 5 ],
	    [ 6, 7, 8, 9 ] );
	my @colours = ( "black", "red", "green", "blue",
			"magenta", "darkgray" );
	for (my $g=0; $g<=$#plot_groups; $g++) {
	    my ($miny, $maxy);
	    my $ptitle = "Baselines: ";
	    for (my $pbl=0; $pbl<=$#{$plot_groups[$g]}; $pbl++) {
		my $tmin = min($WVRPhi[$plot_groups[$g]->[$pbl]]);
		my $tmax = max($WVRPhi[$plot_groups[$g]->[$pbl]]);
		if ($pbl == 0) {
		    $miny = $tmin;
		    $maxy = $tmax;
		} else {
		    $miny = ($tmin < $miny) ? $tmin : $miny;
		    $maxy = ($tmax > $maxy) ? $tmax : $maxy;
		}
		if ($pbl > 0) {
		    $ptitle .= ", ";
		}
		$ptitle .= $baselines[$plot_groups[$g]->[$pbl]]->[0]."".
		    $baselines[$plot_groups[$g]->[$pbl]]->[1]."=".
		    $colours[$pbl];
	    }
	    $offsets_plt->env(min($plt_x), max($plt_x), $miny, $maxy);
	    $offsets_plt->label_axes("Time since midnight UT [m]", "Phase [Deg]",
				     "(".$obsorder[$oo]->{'source'}.") ".
				     $ptitle, { 'colour' => 1,
						'charsize' => 2.5,
						'panel' => ($g + 1) });
	    for (my $pbl=0; $pbl<=$#{$plot_groups[$g]}; $pbl++) {
		$offsets_plt->line($plt_x, $WVRPhi[$plot_groups[$g]->[$pbl]], {
		    'colour' => $colours[$pbl], 'linewidth' => 2,
		    'panel' => ($g + 1) } );
	    }
	}

	# CHECKED HERE, VERY SIMILAR TO wvr_phase.pro

	# Now put out the per antenna phase corrections.
	# Negate them.
	my @p_ant;
	for (my $a=1; $a<=6; $a++) {
	    my $ps = 0.0;
	    my $nant = 6 - ($#excluded_ants + 1);
	    for (my $abl=0; $abl<=$#baselines; $abl++) {
#		my $include_ant = 1;
		if ($baselines[$abl][0] == $a ||
		    $baselines[$abl][1] == $a) {
#		    for (my $cinc=0; $cinc<=$#baseline_exclude; $cinc++) {
#			if ($abl == $baseline_exclude[$cinc]) {
#			    $include_ant = 0;
#			    last;
#			}
#		    }
#		    if (!$include_ant) {
#			next;
#		    }
#		    $nbase++;
		    if (&ant_excluded($baselines[$abl][0]) ||
			&ant_excluded($baselines[$abl][1])) {
			next;
		    }
		    if ($baselines[$abl][0] == $a) {
			$ps -= $WVRPhi[$abl];
		    } else {
			$ps += $WVRPhi[$abl];
		    } 
		}
	    }
	    $p_ant[$a - 1] = -1.0 * $ps / $nant;
	}

	# Write out the corrections to file now.
	my @outstrings;
	for (my $i=0; $i<$wvrTime->nelem; $i++) {
	    my ($day, $mo, $year, $ut) = mjd2cal($wvrTime->index($i) + $df_offset);
	    my ($sign, $hh, $mm, $ss) = time2hms($ut, 'H', 1);
	    $outstrings[$i] = sprintf "%02d%3s%02d:%02d:%02d:%04.1f",
	    $year, $mons[$mo - 1], $day, $hh, $mm, $ss;
	}
	printf OUT "%20s %12.1f %12.1f %12.1f %12.1f %12.1f %12.1f\n", 
	$outstrings[$_], $p_ant[0]->index($_), $p_ant[1]->index($_),
	$p_ant[2]->index($_), $p_ant[3]->index($_), $p_ant[4]->index($_),
	$p_ant[5]->index($_) foreach (0 .. $wvrTime->nelem-1);
	
    }
    print STDERR "done.\n";
    $offsets_plt->close();
    close(OUT);
}

print STDERR "END: wvr_ops.pl\n";

exit;

sub aminmax {
    my $aref = shift;
    my ($min, $max) = $aref->[0];
    for (my $i=1; $i<=$#{$aref}; $i++) {
	$min = ($aref->[$i] < $min) ? $aref->[$i] : $min;
	$max = ($aref->[$i] > $max) ? $aref->[$i] : $max;
    }

    return ($min, $max);
}

sub partial_minmax {
    my $aref = shift;
    my $i = shift;
    my $jlist = shift;
    my $kmin = shift;
    my $kmax = shift;

    my ($min, $max);
    my $z = $aref->[$i];
    $min = $max = $z->[$jlist->[0]]->[$kmin];
    for (my $j=0; $j<=$#{$jlist}; $j++) {
	for (my $k=$kmin; $k<=$kmax; $k++) {
	    $min = ($z->[$jlist->[$j]]->[$k] < $min) ?
		$z->[$jlist->[$j]]->[$k] : $min;
	    $max = ($z->[$jlist->[$j]]->[$k] > $max) ?
		$z->[$jlist->[$j]]->[$k] : $max;
	}
    }

    return ($min, $max);
}

sub check_path {
    my $cmd = shift;

    if (`which $cmd` ne "") {
	return 1;
    } else {
	return 0;
    }
}

sub unique_push {
    my $aref = shift;
    my $nval = shift;

    my $f = 0;
    for (my $i=0; $i<=$#{$aref}; $i++) {
	if ($aref->[$i] eq $nval) {
	    $f = 1;
	    last;
	}
    }

    if ($f == 0) {
	push @{$aref}, $nval;
    }
}

sub mjd2unix {
    my $mjd = shift;

    my $r = bat2unixtime(mjd2bat($mjd));
    return $r;
}

sub which_source {
    # Returns the source being observed at a particular MJD.
    my $mjd = shift;
    my $listref = shift;

    my $mindiff = 2; # in days.
    my $csource = "";

    foreach my $s (keys %{$listref}) {
	for (my $t=0; $t<=$#{$listref->{$s}->{'start_times'}}; $t++) {
	    my $d = $mjd - $listref->{$s}->{'start_times'}->[$t];
	    if ($d > 0 && $d < $mindiff) {
		$mindiff = $d;
		$csource = $s;
	    }
	}
    }
    
    return $csource;
}

sub ant_excluded {
    my $antnum = shift;
    for (my $i=0; $i<=$#excluded_ants; $i++) {
	if ($excluded_ants[$i] == $antnum) {
	    return 1;
	}
    }
    return 0;
}

sub print_help {
    # Print out a help message for the use of this script.
    print STDERR "\n\nPROGRAM: wvr_ops.pl\n\n";
    print STDERR "AUTHOR: Jamie Stevens (Jamie.Stevens\@csiro.au)\n";
    print STDERR "        This program is free software, licensed under GPLv3.\n\n";
    print STDERR " This script is to calculate phase corrections using data from\n";
    print STDERR " the Water Vapour Radiometers (WVRs) installed on the mm receivers\n";
    print STDERR " on the Australia Telescope Compact Array.\n\n";
    print STDERR " Full instructions on how to use the WVR systems and this script\n";
    print STDERR " can be found at:\n  http://www.narrabri.atnf.csiro.au/observing/WVRs.html\n\n";
    print STDERR " The simplest way to use this script is to load your CABB data\n";
    print STDERR " in using Miriad's atlod task, and then run this script with\n";
    print STDERR " the name of that output dataset as the only argument.\n";
    print STDERR "  eg. wvr_ops.pl data.uv\n\n";
    print STDERR " The output of this script is two files for each central frequency\n";
    print STDERR " that it finds in the dataset. One of the files is a PS with\n";
    print STDERR " plots of WVR phase as a function of time and baseline. The other\n";
    print STDERR " is a file with the antenna-based phase corrections that can be\n";
    print STDERR " applied to the data with the Miriad task atwvr.\n\n";
    print STDERR " This script accepts the following options:\n";
    print STDERR "  --cycle\n";
    print STDERR "   -c     By default, the corrections will be output every 5s\n";
    print STDERR "          over the time range of the observation. Use this option\n";
    print STDERR "          to change this cadence. For example, to output the\n";
    print STDERR "          corrections every 10s, use \"--cycle 10\".\n";
    print STDERR "  --refresh\n";
    print STDERR "          When this script runs, it will query the datasets that\n";
    print STDERR "          you specify, and the Narrabri MoniCA server to get the\n";
    print STDERR "          data from the WVRs. Both these operations take a\n";
    print STDERR "          substantial amount of time. This script caches the\n";
    print STDERR "          results of these queries, so that subsequent runs of\n";
    print STDERR "          this script on the same data (to generate the corrections\n";
    print STDERR "          at a difference cadence, for example) will take much\n";
    print STDERR "          less time. To force this script to ignore the cache\n";
    print STDERR "          files and redo the queries, use this option.\n";
    print STDERR "  --wvr-align\n";
    print STDERR "          By default, this script will align the corrections to\n";
    print STDERR "          the times found in the Miriad dataset. To make the\n";
    print STDERR "          script align the corrections to the times found from\n";
    print STDERR "          the WVR data, use this option.\n";
    print STDERR "  --verbose\n";
    print STDERR "   -v     Make this script output more information to the screen.\n";
    print STDERR "  --rpfits\n";
    print STDERR "          By default, this script expects to be given the names\n";
    print STDERR "          of one or more Miriad datasets for which to compute\n";
    print STDERR "          WVR phase corrections. If you would like instead to\n";
    print STDERR "          provide the names of the RPFITS files covering the\n";
    print STDERR "          observing time, use this option.\n";
    print STDERR "  --exclude\n";
    print STDERR "   -x     Exclude an antenna from the computation of the WVR\n";
    print STDERR "          phase corrections. If you want to exclude more than\n";
    print STDERR "          one antenna, use this option more than once. For\n";
    print STDERR "          example, to exclude both antennas 6 and 4 from the\n";
    print STDERR "          computation, give the options \"-x 6 -x 4\" on the\n";
    print STDERR "          command line.\n";
    print STDERR "  --help\n";
    print STDERR "   -h     Print out this help information and exit.\n";

    exit;
}
